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ABSTRACT 

I propose a new method to study computationally difficult problems. I consider a 
new system, larger than the one I want to simulate. The original system is recovered 
by imposing constrains on the large system. I simulate the large system with the 
hard constrains replaced by soft constrains. I illustrate the method in the case of the 
ferromagnetic Ising model and in the case the three dimensional spin-glass model. 
I show that in both models the phases of the soft problem have the same properties 
as the phases of the original model and that the softened model belongs to the same 
universality class as the original one. 

I show that correlation times are much shorter in the larger soft constrained 
system and that it is computationally advantageous to study it instead of the original 
system. 

This method is quite general and can be applied to many other systems. 
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It is well known that there exists a large class of systems which are particularly hard to study 
both analytically and numerically. Such systems are found in condensed matter physics (spin 
glasses, random field models and more generally disordered systems) but also in combinatorial 
optimization problems or in communication theory (error correcting codes). The common future 
of those systems is the existence of a large number of local minima, separated by large barriers. 
The simulation algorithm gets trapped in those minima. 

In this paper I propose a new method, I call soft annealing, in which I enlarge the space on 
which the original system is defined. The original system is recovered by imposing constrains 
on the large system. I propose to simulate the large system with the hard constrains replaced 
by soft constrains. By enlarging the space on which the problem is defined, it is hoped that one 
can get around the barriers and speed up the dynamics of the algorithm, while retaining all the 
essential properties of the system. This is known to be the case for error correcting codespQ. 
I will show that this is also the case for spin glasses and that it is advantageous to study the 
large system. 

I will illustrate the method with the example of spin models in three dimensional cubic 
lattices, but the method is more general. 

Consider a spin model described by the Hamiltonian 

L 

-H° = J x {x,y,z)a(x,y,z)a(x + l,y,z) + J y (x,y, z)a(x,y, z)a(x,y + l,z) 

x,y,z=l 

+J z (x, y, z)a(x, y, z)a(x, y,z + 1) 

where the a's are Ising spins on a cubic lattice of linear size L. 
Now consider the following new Hamiltonian: 

L 

_ H ne W = J x ( x , y, z)a x (x, y, z)a x (x + 1, y, z) + J y (x, y, z)a y (x, y, z)a y (x, y + l,z) 

x,y,z=l 

+J z {x, y, z)a z (x, y, z)a z (x, y,z + 1)+ 

u{a x (x, y, z)a y (x, y, z) + a x (x, y, z)a z (x, y, z) + a y (x, y, z)a z (x, y, z)) 

jjnew con t a i ns three times more spins than H°. To every a(x, y, z) of the original Hamilto- 
nian H° correspond the three spins a x (x, y, z), cr y (x, y, z) and a z (x, y, z). a x (x, y, z) is coupled 
to its near neighboring spins in the direction x, a y (x,y,z) in the direction y and a z (x,y,z) in 
the direction z. The ferromagnetic coupling u couples together the three types of spin on every 
lattice site. For u = H new reduces to 3L 2 decoupled one dimensional chains. For u — > oo, 
a x (x, y, z) = a y (x, y, z) = a z (x, y, z) and H new reduces to H°. I will call the model described by 
H new the soft constrained model and the original model described by H° the hard constrained 
model. 

If H° has a phase transition at (3 = 1/T = (3 C , we expect the phase diagramme of H new 
to be two dimensional. For small u, H new is essentially one dimensional and there is no phase 
transition. For (3 < j3 c there is again no phase transition, irrespectively of the value of u. In 
general, for (3 < f3 c , there will be a line in the (3 = 1/T, v = u(3 plane, separating the high 
temperature from the low temperature phase. We expect the points along this critical line to 
be on the same universality class, i.e. the critical exponents to remain the same, and equal to 
their value in the model described by H° i.e. the u — > oo limit. 

First I verified that all this is true in the case of the ferromagnetic Ising model in 3 di- 
mensions. It is known that in this case f3 c ~ .22 . I simulated H new for (3 = .25 and different 
values of v. I found a ferromagnetic to paramagnetic phase transition for v = v c = .794 ± .005. 
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By finite size scaling I checked that the critical exponents u and the magnetic susceptibility 
exponent are compatible with the best known values for the 3 dimensional Ising model. 

Much more interesting is the case of the spin glass model because it is well known to be 
computationally hard. I simulated the Edwards Anderson model on a cubic lattice with periodic 
boundary conditions. The couplings J are independent random variables taking the values ±1 
with equal probability. This model has been studied a lot (for a review see reference^). It 
undergoes a phase transition for (3 = .90 and the value of the critical exponent v is v = 1.7±.3(3]. 

I studied the soft version of this model for different values of (3, keeping the ratio u = v/f3 
fixed to u = .75 and for the lattice sizes L = 8, 12, 14 and L = 16. For every size I simulated 
1280 realizations of the couplings. For every realization of the couplings I simulated two copies 
a and r and studied the probability distribution of the overlap 

q = (?T3 ) 51 O 7 *^' f ' z)t x (x, y, z) + a y (x, y, z)r y (x, y, z) + a z (x, y, z)t z {x, y, z)) (1) 



I measured the Binder ratio r((3, L) = < q 4 >/(< q 2 >) 2 where, as usually, a is the average 
of a over the coupling samples, and the results are shown in figure 1. The data for the different 
sizes cross at (3 = .96 ± .05. This shows that there is a phase transition for (3 = [3 C = .96 ± .05. 
Figure 2 shows the same data plotted versus the scaling variable s = (f3 — f] c )L 1 / u with v = 1.7, 
the exponent of the Edwards- Anderson modelJS]. We see that the data for the different sizes 
collapse on the same curve except maybe for the L = 8 data which seem to be slightly apart, 
which indicates corrections to scaling for L = 8. The value r c of r at (3 = (3 C is an universal 
quantity, sometimes called the renormalized coupling constant. We find r c = 1.53 ± .03. For 
the original model r c ~ 1.50. So our results are compatible with both the exponent v and the 
renormalized coupling constant r c taking the same value in the soft model and in the original 
Edwards- Anderson model. We conclude that the hypothesis that the two models belong to the 
same universality class is supported by our data. 

Figure 3 shows the probability distribution of the overlaps P(q) averaged over the samples 
for (3 = 1.05 and for the two sizes L = 8 and L = 14. The similarity with the P{q) of the hard 
constrained model^ is striking. P(q) has a peack at a value of q we call qEA- As the volume 
increases, the peack is sharper and qEA decreases, as in the original Edwards Anderson model. I 
measured also the link-link overlap q e of two real replicas a and r, q e = (1/3L 3 ) Pi&jTiTj (the 
sum runs over the nearest neighbor sites of the lattice) and the joined probability distribution, 
P(q, q e ) averaged over the samples. I found that P(q, q e ) is very peacked around a line q e = 
eo + e±q 2 + e2q 4 , with eo ~ .54, e\ ~ .10 and ei ~ .30. The e^'s are to a good approximation L 
independent. I measured the variance w 2 = < (q e — eo — e±q 2 — e2q A ) 2 >■ I found for (3 = 1.05 
^2 = .0004 for L = 8, w 2 = .00015 for L = 12, w 2 = .0001 for L = 14, and w 2 = .00008 for 
L = 16, i.e. when the volume increases P(q,q e ) goes to a delta function P(q,q e ) = P(q)5(q e — 
eo + e\q 2 + e2<? 4 ) This is the case in the infinite range model and in the case of the hard 
constrained model jSJ Ej. We conclude that the hard and soft constrained models belong to 
the same universality class and that their low temperature phases are extremely similar. 

We expect the soft model to be easier to simulate numerically. To verify this one usually 
measures the spin autocorrelation functions, averaged over the samples, 

C h (t) = (l/L 3 )^< CTi (t )a i (i + t) > 

i 

and 

C s (t) = (1/3L 3 )< J2 of(to)of{to + t) + aUt )af(t + t) + af(t )af(t + t) > 

i 

for the hard and soft constrained models respectively, where t is the number of Monte Carlo 
sweeps over the lattice. I compared C h {t) and C s (t) each one measured at the critical value of 
(3, (3^ = .90 for C h {t) and S C = .97 for C s (t). They are plotted in figure 4. 
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There are three curves in figure 4. 

The lower curve (in black) corresponds to C s {t). It is difficult to distinguish between the 
two other curves because they fall on top of each other. One of the two, the curve in blue, 
corresponds to C h (t). We see that C s (t) indeed decays much faster than C h (t). In order to 
make the comparison more quantitative, I plotted also in figure 4 C s (wt) (the curve in red), 
i.e. I rescaled the "time" in C s (t) by an appropriated scaling factor w. We see that, with this 
rescaling, C s (wt) and C h {t) can be made to coincide. 

The scaling factor w was chosen as follows. C s (t) and C h (t) are decreasing functions of t 
and C s (0) = C h (0) = 1. I measure the number of lattice sweeps t h (c,(3) and t s (c,(3) needed 
for C s and C h fall below the value c, i.e. C h (t h (c,(3)) = c and C s (t s (c,(3)) = c. For c = .40 
and L = 12 I found that the time ratio r ? itica \c = .40) = t h (.40, ^)/t s (A0, (i s c ) ~ 28. In figure 
4 I chose w = rf ltlcal (c = .40). We see that with this rescaling of the Monte Carlo time, the 
spin autocorrelation functions coincide for all 1 < t < 10 6 , i.e. one has to perform 28 times 
more lattice sweeps for the hard constrained system in order for the spins to de-correlate the 
same amount as for the soft system. In order to study the size dependence of this difference 
of the correlation times, I measured rf ltlcal (c = .40) for L = 8 and L = 16. I found that 

r critical( c = A ty = for L = g and r critical( c = _ 4 q) = 39 for L = ^ i g ^ritical ( &nd 

therefore the gain in simulation time) strongly increases with lattice size. It is not clear yet 
whether r ^ rittcal obeys some kind of finite size scaling. 

I conclude that, by softening the constrains, the gain in computer time is large, despite the 
fact that the number of spins in H new is larger by a factor of three. This gain icreases with the 
size of the system. 

The method of softening the constrains presented here is obviously not unique. For example 
one could consider 

H' = JijXiJ + u X XijXj tk X k jXi ti 

i,j plaquettes 

The term proportional to u runs over all the plaquettes, i.e. the elementary squares, of the 
lattice. Xi j corresponds to cjjCJj of the original Hamiltonian and the coupling u imposes the 
constrain XijXj^X^^Xi^ = OiOjOjOkQkViQiVi = 1 One of the advantages of the method used 
in this paper is that all the techniques which speed up simulations like multi-spin coding and 
parallel tempering can easily be implemented in the soft model. 

I have illustrated the method with the notoriously difficult example of the three dimensional 
spin glass model. It would be interesting to apply this method also to other difficult optimization 
problems. 

I would like to thank Andrea Montanari for several discussions. 
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Figure 1: Binder ratio r plotted as a function of (5 = 1/T. Black circles represent data for 
linear size L = 16, red triangles for L = 14, green x's for L = 12 and blue stars for L = 8. 
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Figure 2: Binder ratio r plotted versus the scaling variable s (see text). Black circles represent 
data for linear size L = 16, red squares for L = 14, green x's for L = 12 and blue stars for 
L = 8. 
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Figure 3: Average probability distribution P(q) of the overlaps q for (3 = 1.05 and sizes L = 14 
(stars in red) and L = 8 (crosses in black). 




Figure 4: Spin-spin autocorrelations C (t) and C s (t) at (3^ and /?;? (see text) as a function of 
the number of the lattice Monte Carlo sweeps. 



